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We compute the distribution of the number of negative eigenvalues (the index) for an ensemble 
of Gaussian random matrices, by means of the replica method. This calculation has important 
applications in the context of statistical mechanics of disordered systems, where the second derivative 
of the potential energy (the Hessian) is a random matrix whose negative eigenvalues measure the 
degree of instability of the energy surface. An analysis of the probability distribution of the Hessian 
index is therefore relevant for a geometric characterization of the energy landscape in disordered 
systems. The approach we use here is particularly suitable for this purpose, since it addresses the 
problem without any a priori assumption on the random matrix ensemble and can be naturally 
extended to more realistic, non-Gaussian distributions. 

PACS numbers: 75.10.Nr, 61.43. Fs, 64.70.Pf, 61.20.Lc 



I. INTRODUCTION 

The importance of random matrix theory (RMT) can hardly be overstated. Since its initial development by Wigner 
and Dyson to deal with the spectrum of many-body quantum systems it has found applications in areas of 

physics as diverse as disordered systems, chaos, and quantum gravity, to name just a few Most of the time 

RMT has been used as a very powerful tool for the study of the energy-level fluctuations of quantum systems. In this 
case the matrix that RMT is modeling is of course the quantum Hamiltonian of the system. 

However, there is a different context where RMT can be very useful, namely the study of the statistical properties of 
classical disordered systems. By disordered systems we mean not only those cases where quenched disorder is directly 
present in the Hamiltonian, as in spin-glasses, random field models or neural networks, but also systems whose physical 
behaviour at low temperatures is heavily influenced by the self-induced disorder of their typical configurations, as, 
for example, supercooled liquids and structural glasses. In all these systems the properties of the energy landscape, 
or energy surface, are known to be far from trivial. In particular, the presence of many local minima of the potential 
energy is one of the most distinctive features of this class of systems [Ril. An obvious consequence of this fact 
is that the energy surface displays many extensive regions with unstable negative curvature and therefore has very 
non-trivial stability properties [p|JlO|]. In this context a key object becomes the matrix of the second derivatives of the 
Hamiltonian, normally called Hessian, which encodes all the stability attributes of the energy landscape. 

The study of the statistical properties of the Hessian has been an important issue both in the theory of mean-field 
spin-glasses and in liquid theory. In the former case it is often possible to analyze the Hessian in the stationary points 
of the free-energy, having therefore important information on the shape and stability of the thermodynamic states . 
In liquids, on the other hand, the Hessian of the potential energy is the key object in the context of the instantaneous 
normal modes approach where the average spectrum of the Hessian is directly connected to many physical 

observables of the system. In particular, it has been argued that there exists a deep relation between the diffusion 
properties of a liquid and the negative unstable eigenvalues of the average Hessian Q . 

It is evident that in the above context an application of RMT to the study of the statistical properties of the Hessian 
can be potentially very useful. An important remark is the following: the Hessian is a matrix which in general depends 
on the configuration of the system and possibly also on the quenched disorder, when this is present. The basic idea is 
to derive from the distribution of the configurations and from the distribution of the disorder an effective probability 
distribution for the Hessian, which can then be studied in the context of RMT (a recent example of this strategy can 
be found in Hi. 
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Besides, it is clear from the former discussion that an important issue is the analysis of the negative eigenvalues of 
the Hessian, since their presence is related to regions of unstable negative curvature of the energy surface and thus 
possibly to the boundaries of different basins of attractions in the phase space. In particular, the number of negative 
eigenvalues of the Hessian, called the index, is the hrst and easiest measure of instability. As a consequence, all the 
tools devised for the investigation of the index in RMT are particularly relevant in the context of statistical mechanics 
of disordered systems. The average value of the index is trivially related to the average spectrum of the Hessian by a 
simple integration. On the other hand, a more interesting and less trivial quantity is the probability distribution of the 
index. Indeed, while the average index gives a measure of the overall degree of instability of the energy surface, the 
knowledge of the fluctuations of the index around its average value allows a more profound and complete geometric 
description the energy landscape. 

In this paper we compute the probability distribution of the index for an ensemble of Gaussian random matrices 
with a diagonal shift. This ensemble provides the simplest possible model for the Hessian of a disordered system at 
a given energy and represents the ideal context where to develop the technical aspects of this kind of computation. 
Moreover, in the Gaussian context we are able to give non-trivial physical interpretations of our results. 

In order to compute the index distribution we use a fermionic replica method. In the past the replica method 
has been applied to recover standard results in RMT, with variable success. Recently the interest of the community 
has focused again on this method |jlj-|l6| and some indications of the mathematical consistency of the method have 
been provided, even if some strong criticisms still persist |l7|] . The present computation offers an interesting example 
where the replica method can be applied to obtain exact results which are not easily available in the standard RMT 
literature. 

There is also another important reason for using the replica method in the computation of the index, which is 
related to the physical relevance of the Hessian discussed above. As we have seen, RMT can be used once an effective 
probability distribution for the Hessian has been worked out from the distribution of the configurations and from the 
distribution of the quenched disorder. This effective distribution will not be Gaussian in general (unless we consider 
some very particular models) and typically it will not belong to the standard ensembles considered by ordinary RMT. 
By means of the replica method we have in principle no need to assume any specific form of the distribution. 

The paper is organized as follows. In Sec. II we compute the average determinant for matrices of the Gaussian 
Orthogonal Ensemble as a warm-up exercise to fix notation and ideas. We then proceed in Sec. Ill to the main part 
of the paper, where we calculate the average index distribution by means of the replica method, in the limit of large 
matrices. In Sec. IV we apply the previous analysis to the specific case of a mean-field spin-glass model, where the 
Hessian is exactly a Gaussian random matrix. Finally in Sec. V we discuss the general relevance of our results and 
state our conclusions . Technical details of the calculation and the contribution of replica symmetry broken solutions 
are contained in two appendices. 



II. A PRELIMINARY CALCULATION 

Consider the matrix 

M ij = J ij -ES ij , (1) 
where Jjj is an iV-dimensional real and symmetric random matrix with the Gaussian distribution function 

P[J]=2-^ 2 (^) JV2/2 exp(-^TrJ^. (2) 

We have introduced a diagonal shift E in order to mimic what in general happens in disordered systems, where M 
represents the Hessian of the Hamiltonian. In this context we expect to find very few negative eigenvalues of M at 
low energies, because of the dominance of minima at very low energies. This is the effect of the shift in (jlj) and we 
therefore shall refer in the following to E as to the energy. 

The average density of eigenvalues, or spectrum, of M is defined by 

p(X; E) = --L ImTr (A - M + i e)" 1 = Im JL logdet(A - M + ie) , (3) 
7r7V irN oX 

where the bar indicates the average over distribution (g). It is well known that for the Gaussian ensembles the 
spectrum p in the limit N — > oo is given by a semi-circle centered around A = —E, that is 
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p{\- 1 E) = -^A-{\ + EY , (4) 

Z7T 

while p is zero outside the semi-circle support ||. 

In order to fix our notation and to acquire some familiarity with the method we will use, we compute in this section 
the average determinant of M . In general this is not a self-averaging quantity, in the sense that fluctuations around 
the mean value do not decrease in the limit N — > oo. The correct object to average is in principle the logarithm of 
the determinant, as it appears in the definition of p, since this is an extensive quantity. However, it is a particular 
property of the Gaussian case that the determinant is self-averaging at the leading order, so that the calculation of 
det M is an interesting and simple warm-up exercise for what we want to show later. 

We can write the determinant by means of a Gaussian integral over TV-dimensional fermionic vectors (ip, ip) 



N 



(5) 



det M = J dip dip exp 
We now average over the symmetric matrix Jki 

/ N N \ 

(6) 

To decouple the quartic term in the fermions we perform a Hubbard-Stratonovich transformation 

(7) 



r ( N - 1 " _ _ \ 

det Af = / dtp dip exp E 2j^M>j _ ojV X/ ^i^i^i 

J \ i=l i,j=l J 



r _ ( N - N - \ 

det M = j dip dip dq exp I E ip{pi — — q 2 + iq ^ ip^rpi J 



and after integrating out the fermions we obtain 

detl7= / dqe NS{q > E) , (8) 



with 
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S(q,E) = -- q 2 +log(-E-iq) . (9) 

This integral can be solved exactly in the limit TV — > oo by means of the steepest descent method. The procedure is 
quite standard jL8|, but we briefly summarize it for the sake of clarity. In order to calculate integral (|J) in the large 
N limit we must select a path of integration 7 in the complex plane, which satisfies the following conditions: 

i) The integral along 7 must be equal to the integral along the original integration path (in our case the real axis). 

ii) The imaginary part of the action S(z, E) (or phase) must be constant along 7. 

Hi) The path 7 must pass through at least one of the saddle points of the action S(z, E). 

The integral along 7 can then be computed using the Laplace method and it is given, at the leading order, by 
the integrand evaluated in the maximum of the real part of S along 7, that is, in the saddle point of the whole action. 
In the case where many maxima lie on 7, only those with the largest real part of S contribute to the total integral. 

In our case the action S has two saddle points in the complex plane, given by 

q ± = ^E±^4-E 2 . (10) 
The regions of constant phase passing through q+ and q_ arc defined by 

7+ : ImS(z) = lmS(q+) 

7_ : Im5(z) = ImS(q-) . (11) 

These regions satisfy by definition conditions (ii) and (Hi) and thus the correct path of integration 7 must be built 
by using the different branches of 7+ and 7_ in such a way to satisfy condition (i). We can distinguish three different 
regimes: 
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• E < —2: For these values of the energy the imaginary part of the action is the same for the two saddle points. 
The constant phase region is shown in Fig.l: it is clear that there is only one path 7 satisfying condition (i) which 
can be built by means of the different branches of the constant phase region. This path is almost parallel to the 
real axis and passes through q + , but not through Indeed, the path parallel to the imaginary axis, which passes 
through both the stationary points, does not conserve the original integral. The only stationary point contributing 
to the integral is therefore q + and we have 

tetM = e NS( q+ ,E) =2 -N ^eT—^ N e N(\E\ + VE^I) 2 /8 (12) 

In this energy regime the spectrum p has support completely contained in the positive semi-axis and we thus expect 
the average determinant to be positive, as it is. 



Constant phase region — 
Real axis 



FIG. 1. E < —2: The region of constant phase (dashed line) and the real axis (full line). The two small circles indicate the 
positions of the two saddle points, q + (up) and g_ (down). The only suitable path of integration 7 passes just through q + , 
since the original integral is not conserved on the orthogonal path. The case E > +2 is specular to this one. 



• E > 2: The support of the spectrum is now entirely contained in the negative semi-axis, so we expect all 
eigenvalues of the matrix to be negative. In this case the path 7 passes only through the saddle point <?_, and we 
thus find for the determinant 



ditM = e JVS («- £ ) = (-l) JV 2- JV (\E\-y/E^Z) e N(\E\+VE^if/^ (13) 

with the correct prefactor (— l) N indicating that all eigenvalues are negative. 

• — 2 < E < +2: In this regime the situation is very different. In Fig. 2 we plot the region of constant phase: the 
only path 7 which satisfies condition (i), passes now through both the saddle points q + and </_. It must be noted that 
in this case the imaginary part of S is different in q + and g_, so that actually the global region of constant phase 
plotted in Fig. 2 is the union of two different regions, 7+ and 7_. On the other hand, the real part of S is the same 
in the two stationary points, and therefore they both contribute to the integral. We have 

fc^M = e NS(q+,E) + e NS(q-,E) = ^_^N a (E) e N(E 2 -2) /i+log 2 ^ ^ 

where 




(E) = -arctg ^— + —E\J A - E 2 , (15) 
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Constant phase region 

Real axis 




FIG. 2. —2 < E < 2: The region of constant phase (dashed line) and the real axis (full line). The two small circles indicate 
the positions of the two saddle points, q+ (right) and q_ (left). In this case the correct path of integration 7 passes through 
both the saddle points. 



At these values of the energy the spectrum of M is partly contained in the negative semi-axis, so that a non-trivial 
fraction of the eigenvalues is negative. The interesting point is that the interplay between the two saddle points gives 
rise to the correct sign of the determinant. Indeed, it is easy to check that a(E) is exactly the mean fraction of 
negative eigenvalues of M, that is 

a{E) = ( dXp{X;E) . (16) 



Note that the mechanism we have described above, given by the interplay between the two saddle points and the 
paths of integration, is crucial in order to obtain the correct result for the determinant of M. 

III. THE INDEX DISTRIBUTION 

The index Xm of a matrix M , defined as the number of its negative eigenvalues, can be computed from the following 
formula (n]] 

T M = — lim [logdet(M - ie) - logdetfM + ie)] . (17) 

Z7TZ e^O 

The meaning of this relation is quite clear: the function f(z) = logdet(Af — z) has a cut on the real axis at each 
eigenvalue of M, such that by means of the limit in ( |l7j) we are crossing as many cuts as negative eigenvalues are 
present. Besides, this formula can be simply obtained by integrating the non-averaged spectrum (|^) from minus 
infinity up to zero. In the case we are considering, the index is a function of the energy E and its average value is 
given by Na(E) (Eq.(|l|)). 

We are interested in calculating the average probability distribution of the index, at a given energy E, that is the 
probability P{K] E) to have a matrix M with index Tm equal to K, at energy E 



P(K;E) = S(K-X M (E)) . (18) 

In the following it will be important to distinguish between the extensive index K , which is a positive integer between 
and N, and the intensive one k — K/N, which takes values in the continuous interval [0, 1], and whose probability 
distribution is 
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p{k;E) = 5(k-l M (E)/N) = NP(Nk; E) . (19) 



Note that the limit N — ► oo is well defined only for p(k; E). 
From Eqs.Q and © we S et 
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P(Jf; E) = ^J dp e-* K G(ji, E) , (20) 
where 

G(m, £) = det" /2,r (M - »e) det^ /27r (M + ie) . (21) 

We now make use of the replica method to represent the powers of the determinants in G(/x, E) as analytic continua- 
tions of integer powers 

det ±M/2w (Af T ie) = lim det n± (M=Fie)- (22) 
n± — >±/i/2-7r 

By introducing two different sets of iV-dimensional fermionic vectors (x±, X±) with r = 1, . . . , n± , we can rewrite the 
determinants as 

det^/^M T ie) = lim , / ^xi^xi exp f - |S±(M T ie)*£ J , (23) 

where the sums over the matrix indices i,j are hereafter always understood. We can write everything in a more 
compact fashion by introducing the Grassmann vectors (f/> , ip a ), with a = 1, . . . , (n+ + n_), defined as (see also [p0|) 

(tpi, . . ■>(„++„_)) = (x+, ■ ■ -,X+ + ,X-, ■■■,X-~) , (24) 

together with the matrix 

e ab = diag(e, . ..,£, -6, . .. ,-e ) . (25) 
Note that both ?/> a and e a (, have replica dimension n = {n + + rt_) — ► 0. In this way we have for G 



G(fi,E)= lim / DV-DV'exp 

Tl±— >±/n/27T J 



afc - i£ab) V'b 



x&=l 



(26) 



The average of G over the distribution of J can be computed by a generalization of the procedure of the previous 
section, the main difference being the fact that we have an extra replica dimension, so that the variable q must be 
replaced by a matrix Q a b- For the sake of completeness the details of the computation are in Appendix A. We obtain 



GfaE) = J DQ 



e NS(Q,E) f (2?) 



with 



S(Q, E) = - ilr Q 2 + log det(-25 - iQ) , (28) 

and E a b = E5 a b + itab- Note the similarity with equations (||) and (|^). The matrix Q is an n x n self-dual real- 
quaternion matrix |2l],[l5| (see Appendix A). It has 2n 2 — n degrees of freedom, and is diagonalized by transformations 
of the simplectic group Sp(n). In ( p8f ) we see for the first time the role of e as a symmetry breaking field. The matrix 
E has an upper block of size n + which contains +ie and a lower one of size n_ with — ie, so that the action is only 
invariant under Sp(n) / ' Sp(n+) x Sp(n-), and the full invariance under Sp(n) is only recovered in the limit e — > 0. 
However, how exactly the symmetry breaking affects the calculation will become clearer below. 

We can evaluate the integral (|27j) by means of the steepest descent, or saddle point, method, which becomes exact 
for large N. The saddle point equation for the matrix Q reads 
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Q = iiE + iQ)- 1 



This equation can be solved assuming for Q a diagonal form, Q a b — z a b a b- We have two different sets of equations, 
one set for the elements belonging to the upper block, zi u \ and a second set for the elements of the lower block, z$p . 
The only difference between the two sets is, of course, the sign of e, 

4 M) =i (e + ie+ iz { a u ^ upper block 

2« =i (E-ie+ i4°) 1 lower block (29) 



Each one of these two sets of equations has two solutions, z^P for the upper block, z± for the lower one, namely 



(«) 

4 

.(0 



-(E + ie)±-y/4-(E + ie)* , 
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4> = -(E-ie)±- y/4-(E-ie)* . (30) 
For all values of the energy such that |4 — E 2 \ 3> e these solutions can be expanded in powers of e and read 

'±=^ + f % ±l UT=w) +0{c2) - (31) 

where q± are given in equation (|To|). 

There are some important things to note here, related to the fact that the presence of e crucially modifies the 
mutual relevance of the different saddle points. We have seen in the previous section that in the regime —2<E<2 
the correct integration path 7 passes through both the saddle points q+ and g_ (see Fig. 2). This is true also in the 
present case, when a value e ^ is considered: for each z a the path 7 passes through z+ and z- and in principle 
both the saddle points must be taken into account. However, when we look at the real part of the action S, we now 
discover that the contribution of one saddle point is exponentially dominant over the other by a factor exp(— Ne). 
This is in contrast with the case of the previous section, where the real part of S was the same in the two saddle 
points. 

The crucial point is that, due to opposite sign of e in the upper and lower blocks, the real part of the action is 
tilted in opposite ways in the two blocks and, as a consequence, the dominant saddle point becomes z + for the upper 
block and z_ for the lower one. We now start to understand the way in which e works as a symmetry breaking field: 
without e the two saddle points have the same weight in the integral and we have to consider both of them. With 
e, the weights are modified in opposite ways for the upper and lower blocks. In order to apply the steepest descent 
method we must perform the limit N — * 00 before the limit e — > 0, and this selects just one different saddle point for 
each of the two different blocks, dumping completely the non-dominant contribution. As a result, when at the end 
e^Owe have selected q+ for the upper block and q- for the lower one. This is very reminiscent of what happens in 
statistical physics, where, in order to break a symmetry by means of an external field, the thermodynamic limit must 
be performed before sending the field to zero. 

On the other hand, for energies \E\ > 2, the effect of e is harmless, there is no qualitative change from the situation 
described in the previous section and the same kind of saddle point for the upper and lower block contributes to the 
integral. 

We can now proceed in our computation. We will focus first on the region — 2 < E < 2, where the typical spectrum 
is not positive defined and where we thus expect a more interesting index distribution. According to the above 
discussion on the dominant saddle points, we must consider the following form for Qsp- 

Q SP = diag( 4 u) , zfi z®,...,z® ) . (32) 

n+ n_ 

This form is invariant under the unbroken group S'p(n+) x Sp(n-) of replica symmetry transformations, and in this 
sense we shall refer to it as a replica symmetric (RS) saddle point |l(],[l4|]. We note that Eq. (pl|) is invariant under 
the simultaneous action of complex conjugation and inversion of /j, which after replicating becomes n± — > n T , and 
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that our saddle point satisfies this invariance. If we plug expression (32) into Eq.(|27j), we obtain after taking the limit 
e -> 



GQi,E) 



-E-iq+f'** {-E-tq_y N ^ e-N^l-il.)/^ 



exp [i/iNa(E)} 



(33) 



where a(E) is the average fraction of negative eigenvalues given by Eq. dig) . From ( P0|) we finally get the probability 
p(k,E) in the limit N — > oo, 



p(fc,£) = «y [jfc - a(E)] 



(34) 



This result is very reasonable, but also rather trivial: the probability distribution of the intensive index is a <5- function 
peaked on its average value in the limit N — > oo. In order to observe a non-trivial behaviour we need to consider the 
scaling with N, that is, the distribution of the index for large but finite N. This is particularly important if we are 
interested in the distribution of the extensive index, as for example in the case of disordered systems, where we want 
to know the change in the probability of different stationary points when variations of the index of order one, not of 
order N, are considered. 

To go beyond result (34), we must consider fluctuations around the saddle point ([52]). The general procedure is 
discussed in Appendix B. As expected there are three kinds of fluctuations: within the upper block, within the lower 
block, and those which mix the two blocks. Their corresponding eigenvalues and degeneracies are, 



i + 4 u) 4 u) =(i+^)- 



i , (0 (0 
1 + z w z y 



(1 + q 2 -) 



= 1 



\A- s2 /4 

y/l-E*/4 
+ Q(e 2 ) 



0(6 2 ) 

+ 0(e 2 ) 



d u = 2n + — 
d\ = 2n 2 _ — 7i- 
dm = in+ri— . 



(35) 



The first two sets of eigenmodes are massive modes, in the sense that their eigenvalues are O(l). The third set are 
soft modes: for vanishing e they would correspond to zero modes associated to the restoration of the Sp(n+ + n_) 
symmetry; for small non-zero e they become soft vibrations. Integrating over the fluctuations, we obtain 



(n%-n+/2) -(n 2 _-n-/2) - 2 n+n 



G(fi, E) = uj u ll-j 
In the replica limit n± — > ±fj,/2n this quantity becomes 



G(/z, E) = exp 



ifiNa(E)- f^log 



fi 2 f V^u^l 



exp [ifj,Na(E)] 



f log p 

47T V Wl 



From Eq. (p0|) we obtain the distribution for the extensive and intensive index for finite but large TV: 

1 ; [K - Na(E) + f3(E)} 2 \ 



P(K,E) 
p(k,E) 



2ttA(E) 



exp 



N 2 



■ exp 



2A(E) J 
N 2 [k~a(E)+(5{E)/N] 



2ttA(E) ""^ \ 2A(E) I 

These are Gaussian distributions peaked on the average value a(E). Indeed the shift, 



0(E) = ^arctg 



E 



.V4-£ 2 , 

is of order one and is not relevant at large enough values of N. The variance A(E) is given by 



A(E) = \ lo, 



that is 



A(J0 = -a log [27r 2 e-Vo(^) 2 ] 



(36) 
(37) 

(38) 
(39) 

(40) 

(41) 
(42) 
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where we have defined po(E) = p(X = 0; E) (see Eq.(^)). This result for the variance can also be obtained by the 
method of orthogonal polynomials where e plays the role of a high frequency cutoff p2| , p3p |] . 

The fact that expression ( ff2f ) still depends on e can seem rather unphysical, especially when we consider the fact 
that the limit e — > has to be performed. However, we have to remember that we are looking at finite N corrections, 
and this very fact makes the parameters e and N no longer independent. In this way the presence of e translates in a 
more physical N dependence and this allows us to compute the scaling of the index distribution with the matrix size 
N. Before discussing the result we have obtained for the index distribution, we have therefore to address the problem 
of the relation between e and N . 

There are mainly two different reasons why e and N are related. First, as we have previously noted, there is a 
precise interplay between the two limits, N — > oo and e — > 0, when the saddle point approximation is used in order 
to solve integral (p7|): the symmetry breaking due to e works only if e — > after N — > oo, as in any thermodynamic 
calculation. If N is kept finite, we need a value of e big enough to guarantee the dominance of one saddle point over the 
other. We have seen that the role of e is to modify the real part of the action in such a way that along the integration 
path 7 one saddle point is weighted more than the other. However, if e is too small, also the non-dominant saddle 
point may give a non-negligible contribution to the integral. To avoid this fact we need the secondary contribution to 
be suppressed also at finite N and to vanish when the limit N — > oo is considered. The suppression factor is given, 
at order e, by 

e -N[S(z™)-S(z™)] = e -2nNep (E) ^ 

for the upper block (for the lower block an analogous expression is valid) . In order for the suppression factor to vanish 
it must hold 

eN -> oo , N -> oo . (44) 

This imposes a lower bound for e when N is finite. A natural general choice is therefore to assume 

e = Ni=m ' (45) 

where the exponent S(N) has to satisfy the relation S(N) log A — > oo. The simplest possibility is, of course, a constant 
value of 5. However, as we shall argue immediately below, this would not be consistent with the second condition we 
have to impose on e. 

The second bound for e comes from the following observation. When we perform our calculation with a finite value 
of A and of e, there are of course two different kinds of corrections to the asymptotic exact result: the first kind is 
related to the saddle point approximation and brings corrections which scale as inverse powers of A. The second is 
related to the non-zero value of e and brings corrections which scales as powers of e. Consistency requires that in the 
final result the error introduced by considering a finite value of e must be of the same order as the terms we discard 
in the expansion in 1/A. It can be easily shown that the corrections to the index distribution (Bq) for finite e are of 
order e 2 , that is 

G{p) = exp(Na(E) + p(E) + 0(e 2 )) . (46) 

On the other hand, by considering the Gaussian fluctuations around the saddle point, we are discarding terms of 
order 1/iV 2 in the exponent of (]46]). Thus, we must impose the condition 

Equation ( f47|) is consistent with equations (ff5|) and (|^|) only if, 

5(N)->0, 6{N)logN-+ oo , N -> oo . (48) 
In this way we finally get for the variance the result, 

A(E, N) = ^ log Utt 2 N^- s ^ Pq (E) 2 } = 1 ~ y } log TV + A log(27rp ) , (49) 

where we have taken e = \/2N 1 ~ 5 , the factor 1/2 being consistent with equation (H) at E = 0. This result agrees 
very well with numerical simulations: in Fig. 3 we plot the variance A as a function of log AT, obtained by exact 
numerical diagonalization. A linear fit gives 
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A = -^-logN + \ log(27rpo) , a = 1.005 ± 0.006 , b = 1.993 ± 0.003 . (50) 

7T 7T 

This same scaling for the variance has been found also in p^ , where a completely different method based on the 
invariance properties of the Gaussian Orthogonal Ensemble and the dominance of intrinsic binary correlations was 
used. In Appendix B we show in details that the contributions of the other possible saddle point solutions of the 
whole integral (|27]) to the index distribution are smaller by inverse powers of log N in this energy region, therefore 
the scaling with N is correctly reproduced by equations (|3|) , ( |39"| ) and Q4S[ ) . 



0.9 




FIG. 3. The variance A as a function of log N for E — 0, obtained by means of exact numerical diagonalization on the 
Gaussian Orthogonal Ensemble. The full line is the linear fit. 



We can finally analyze the significance of our result, equations (|38[) , (|39|) and (|49|), in the energy regime — 2 < E < 2. 
What we see is that the variance of the probability distribution of the intensive index goes to zero in the limit N — oo 
and this was quite expected, given our former result (|34|). On the other hand, the variance of the distribution of the 
extensive index diverges logarithmically for TV — > oo. The meaning of this result is the following: on the one hand the 
probability of finding a matrix with an index density different from the average one, that is with an extensive index 
X Na + O(N), is zero in the limit N — > oo. But, on the other hand, the probability of having a matrix whose 
index differs from the average one for a number of negative eigenvalues of order one, i.e. X ~ Na + 0(1), is exactly 
the same as the probability of having a matrix with the average index, in the limit ./V — > oo. As we shall see in the 
next section, this fact has some very interesting physical consequences in the context of disordered systems. 

Let us now look at the other energy regions. First of all we note that the derivation of equations (|3l|), (|35| ) and 
( [43| ) holds as long as the energy is such that po(E) is of 0(1). But this condition breaks down when the energy gets 
close to ±2 and po(E) « 1. In this region the procedure previously adopted to compute the index distribution has 
to be modified. Indeed, when po(E) becomes too small the suppression mechanism ([h|) starts being inefficient, and 
the saddle point ( |32| ) is no longer the only one contributing to the integral. At some point the excitations which were 
treated as soft modes in (|35| ) must be considered as zero modes connecting equivalent saddle points: there exists a 
manifold Sp(n + + n_) / Sp(n+) x Sp(n^) of saddle points and the original replica symmetry under Sp(n + + n_) is 
restored. At this stage e plays no longer any role and it can be taken to zero. The massive modes are the same as in 
Eq. d35| ) , and after integrating over them and exactly over the degrees of freedom associated with the zero modes we 
obtain (up to trivial factors in the replica limit) 



G{i^E) = N^-V n n + ++n _ a," (n +- +/2) cr (n -" n - /2) exp [i^Na(E)] , (51) 

where V™^ + „_ corresponds to the volume of the manifold of saddle points solutions (see Appendix B). At this point one 
has to analytically continue the previous expression for ?i+ — * (i/2ir, n — > 0. The volume V™^ +n _ is finite for rt+ = 
and it is zero for positive integers [Q. Its analytic continuation is an oscillatory function of n+, with exponentially 
increasing amplitude [T^l , so that the presence of this factor in the former equation makes the index distribution 
non-Gaussian. However, as long as this analytic continuation is finite for non-integer n + ~ 1, the distribution can be 
approximated by a Gaussian with variance 
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A(E ~ ±2) = — log [Sir 4 N po(E) 3 ] 



(52) 



We can see from ( p2[ ) that the variance still scales as log TV. However, when \E — 2| ~ l/iVi , we have po(E) ~ I/AT 3 
and a further crossover takes place: the variance A(E) becomes of order one meaning that the index distribution is 
dramatically more peaked around its typical value as we approach E = ±2. Note also that when E ~ — 2 + l/N? 
the typical index a(E) becomes of order 1, meaning that in this region matrices with O(l) negative eigenvalues are 
dominant. Summarizing, in the energy regime where the average number of negative eigenvalues is of order one, the 
fluctuations around the mean value become of order one too. 

When the energy is exactly at the threshold values E = ±2 we have a special case since the saddle point equations 
for the eigenvalues have a single degenerate solution, and the harmonic terms in the expansion around the saddle 
point vanish. It is not difficult to show that the distributions here become 

P(K, E -> -2+) = A" 1 8{K) , P(K,E->2~) = N- 1 S(K - N) . (53) 

The calculation in the regions \E\ > 2 is completely straightforward since e plays no role from the beginning. As 
mentioned before, the same kind of saddle point has to go in both blocks, so that we have 

g SP =diag(zi" ) ,...,zi u) ,zi i) ,...,zi i) ), (54) 

n+ n_ 

where the plus (minus) sign corresponds to negative (positive) energies. There is only one kind of massive fluctuation 
with degeneracy 2n 2 — n, which goes to zero in the replica limit, and thus the integration over fluctuations gives a 
trivial prefactor. The final result for the distribution of K is 

P(K E)-S N ~ l5 W E< ~ 2 (55 ) 

which coincides with the limiting behaviour ( |53"| ) of the distribution in the region —2 < E < 2. Thus, while in the 
energy region —2 < E < 2 values of the index with an O(l) difference from the typical one have a finite probability, 
here the index distribution is so much peaked on the typical value that even small changes in the index have zero 
probability. 



IV. AN APPLICATION TO DISORDERED SYSTEMS 



In this section we consider a mean-field spin-glass model, that has been extensively studied in the last years and 
whose thermodynamical as well as dynamical features are very well known, namely the p-spin spherical model [pi|-p8| . 
Our aim is to use the results of the calculation we have carried out in the previous section, in order to have a better 
understanding of the statistical and geometrical properties of the energy landscape for this model. 

This problem is by itself relevant, because both the static properties and the peculiar off-equilibrium dynamical 
behaviour of mean field spin-glasses, and in particular of this model, are known to be deeply related to the distribution 
of the minima and of the saddles of the Hamiltonian ||7|,p| jic|j20tl . Moreover, it is now commonly accepted that the 
p-spin spherical model shares many common features with structural glasses, which are presently one of the major 
challenges for statistical mechanics. Indeed, notwithstanding the completely different form of the Hamiltonians, some 
structural glasses (in particular fragile glasses) and the p-spin spherical model have a very similar structure of the 
energy landscape p9j| . Therefore, a thorough investigation of the energy landscape for the p-spin spherical model is 
important also for a better understanding of structural glasses. 

As already stated in the Introduction, knowing the index distribution of the Hessian at various energies is equivalent 
to knowing the fluctuations in the stability of the energy surface. In other words, the index distribution tells us what 
are the dominant stationary points of the Hamiltonian (or saddles) at a given energy, and, more importantly, what is 
the probability distribution around the typical saddles, thus providing an insight on the mutual entropic accessibility 
of different stationary points. This is what we are going to describe in this last section. 

The reason why the p-spin spherical model is particularly appropriate for an application of the above calculation 
and concepts is the following: when we look at the stationary points of the Hamiltonian of this system, we find that 
the Hessian matrix M in such stationary points, behaves as a Gaussian random matrix of the same kind as the ones 
considered in the calculations above. More specifically, if we classify the stationary points of the Hamiltonian in terms 
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of their energy density E, we find that the Hessian M(E) in these stationary points is a random matrix of the form 
(see for instance p7|) 



Mij (E) = Jij-ESi 



(56) 



where Jy is an iV-dimensional real and symmetric random matrix with the same Gaussian distribution as (g|), and 
where N is the size of the system. The spectrum of the Hessian in the stationary points is therefore, 



where E t h is the so-called threshold energy, which depends on the parameters of the model (in the previous sections 
it was \E t h\ = 2). Given the particular shape of the Hessian, we can completely disregard the details of the p-spin 
spherical model and assume the results obtained in our calculation as the starting point, interpreting these results in 
terms of probability distributions of the stationary points of the Hamiltonian. 

Let us begin our geometric analysis of the energy landscape from very low energies. When E < — \Eth\ the semi- 
circle is entirely contained in the positive semi-axis and the average determinant of the Hessian is positive: this is 
the region dominated by minima, as the index distribution ( p5| ) shows. Moreover, as we have already noted in the 
previous section, the probability of finding a stationary point with an index different from the typical one (i.e. 0) 
is zero. Minima are strongly dominant in this energy regime. A more careful analysis p0| shows that even in this 
regime there are saddles with non-zero index, but the probability of these objects is exponentially small in N, that is 



This result is obtained by considering non-symmetric contribution to the saddle point equations (see ]10| ] and Appendix 
B) and, consistently with equation (p5[), it gives a contribution too small to be caught by simply analyzing fluctuations 
around the dominant saddle point. The above result shows that at low energies minima are exponentially dominant 
over saddles of order one, and even more dominant over saddles with extensive index. In this sense we shall call 
this region the decoupling regime, since at any energy below the threshold only one kind of stationary points, namely 
minima, dominates. When we raise the energy, we finally arrive at E = — \E t h\'- here the semi-circle touches the 
zero and the decoupling between different stationary points is no longer true. Indeed, it can be proved []l0| that 
Q{Eth) — 0, meaning that at the threshold energy minima and saddles of order one have the same probability. 

Thanks to the calculation of the previous section we are now in the position to answer the following question: 
What happens above the threshold energy ? From a simple inspection of the semi-circle law it is clear that above 
the threshold saddles become important, since many negative eigenvalues appear and the average index Na(E) is 
non-zero. Yet, in order to have information on the degree of decoupling of the stationary points, the simple typical 
index Na(E) is not enough. The reason is the following: the knowledge of the typical index does not tell us whether 
at that same energy other stationary points, different from the typical ones, do or do not have non-zero probability. In 
this sense the mutual entropic accessibility of different stationary points is encoded in the index distribution P(K, E), 
which reveals to what extent the typical saddles are dominant over the non-typical ones. 

From equation ( |38| ) we see that in this regime not only the dominant stationary points are saddles of order N, but, 
also, that the probability of finding a minimum is of order e~ N . The decoupling between minima and dominant 
saddles is therefore much more dramatic than the one we found below the threshold. On the other hand, because 
of the divergence of the variance A with N (equation (fl9])), we see that there is a mixing among saddles with the 
same intensive index: the probability of having a saddle whose index differs from the average by a number of order 
one, is the same as the probability of the typical saddles |30|. In other words, the main result is that there is no 
decoupling among saddles with the same intensive index, so that a mixing of different stationary points occurs, while 
still a decoupling exists between dominant saddles and minima. 

Summarizing, we can therefore distinguish two energy regimes where the probability distribution of the stationary 
points, and therefore the geometric structure of the energy landscape, is very different: a decoupled regime for E < E t h 
and a mixed regime for E > E t h- Interestingly enough, the threshold energy E t h is exactly the asymptotic energy 
where a purely dynamical transition occurs: below a critical temperature Td, the ergodicity is broken and the system 
is no longer able to visit the entire phase space in its time evolution, remaining confined to an energy level higher 
than the equilibrium one. This 'dynamical energy' is equal to E t h [p6|-p8[. 

This suggests us to relate the information we have on the distribution of the stationary points, following from 
the index distribution, to the dynamical physical behaviour of the system. Above Td the equilibrium energy E of 
the system is higher than threshold value E t h and therefore belongs to the mixed regime: the equilibrium landscape 




(57) 



P(K, E) 



e -KNSl(E) 



K = 1,2,3,... . 



(58) 
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explored by the system is dominated by saddles of order N which, as we have shown, are all equally relevant up to 
variations of the index of order one. This means that all these unstable stationary points are equally accessible to 
the system in its time evolution. As Td is approached the equilibrium energy E gets closer and closer to E t h, and 
the properties of the equilibrium landscape change accordingly to the behaviour of P(K, E) we have discussed in 
the previous section: when E ~ — \Eth\ + 1/^V 2 / 3 saddles with index of order one become the most relevant and the 
variance of the index distribution is now finite. This means that minima start having a finite probability in this energy 
regime. The range of temperatures where this behaviour takes places is of order 1/iV 2 / 3 and shrinks to zero in the 
thermodynamic limit. Below Td, the equilibrium energy belongs to the decoupled regime, that is E < Eth' minima 
are now dominant and saddles of any order have exponentially vanishing probability. We can therefore interpret Td 
as the temperature where a geometric transition occurs from a regime of strong mixing of the stationary points to a 
regime of equally strong decoupling. 

V. CONCLUSIONS 

In this paper we computed the average index distribution for an ensemble of Gaussian random matrices. We find a 
result which is in optimum agreement with exact numerical diagonalization. This computation is, in our opinion, an 
interesting example where the fermionic replica method, together with a careful asymptotic expansion of the integrals, 
gives correct results. We hope that the present work can therefore contribute to clarify the role of the replica method 
in the context of RMT. 

Besides, and this was our main purpose, the index distribution provides a really useful tool for investigating the 
geometric structure of the energy landscape in disordered systems. In the previous section we applied this tool to 
the simple case of the p-spin spherical model and discussed the physical consequences of our results. In general, 
the task of computing the distribution of the index of the Hessian is not as simple as in the p-spin model. The 
main reason is that the Hessian usually does not behave as a Gaussian random matrix, because, as noted in the 
Introduction, its distribution is determined both by the distribution of the quenched disorder and by the distribution 
of the configurations. However, the same procedure we adopted in this paper can also be applied to these more 
complicated cases, with the appropriate modifications: to compute the index distribution at a given energy E, one 
has to average over the distribution of the disorder and integrate over the relevant configurations belonging to the 
manifold of energy E Jh}. This is the reason why the method presented in this paper is particularly suitable for this 
task, since it addresses the problem without assuming any particular form for the distribution of the Hessian. 

Finally, there have been recently some attempts to find a connection between the occurrence of a thermodynamical 
phase transition and the change in the topology of the configuration space visited by the system at equilibrium |nj . 
For various non-disordered models which present a second order phase transition it has been shown via numerical 
simulations that the fluctuations of the curvature of the configuration space exhibit a singular behaviour at the 
transition point. This is similar to the behaviour described in the previous section for the p-spin spherical model, 
where the average fluctuations of the index (fL^ ) at the equilibrium energy encounter a dramatic change as the 
dynamical transition is approached (Tcfl . This suggests first of all that also in disordered systems a connection 
between thermodynamical behaviour and topology of the configuration space exists. Besides, the case of the p-spin 
is also peculiar in this sense: it presents a static phase transition which is thermodynamically of second order, but 
it is discontinuous in the order parameter pi| and exhibits a purely dynamical transition at a higher temperature 
p6f . As we have shown, in this case a dramatic change of geometrical properties occurs at the dynamical transition, 
indicating that a more complex situation probably holds for disordered systems which present this sort of behaviour. 
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APPENDIX A: 



In this Appendix we give for completeness the standard procedure used to average and integrate out the fermions 
which gives the sigma model (p7|)-(p8|) |2lf . 

The average over the Gaussian Orthogonal Ensemble (GOE) of Eq.(|26|) yields, 

G{n,E)= lim / DtpD^exp ( E ab ip a ■ ip b + ^tp a ■ ip b ipb ■ 4>a ~ 7^7 V> a ' i>b^ b • i/> a ) , (Al) 

n±->±/i/2-7T J \ ZN ZN J 

where summation over repeated replica indices is implicit and the dot stands for contraction of spatial indices. We 
can define the following n x n matrix whose components are quaternions, 

3 

A ab = A° ab l+J2K^s, (A2) 

8=1 

where 

A ab = ^ a -^ b +^ b ^ a ), (A3) 

Kb = \($ a -i>b-i ) a- i> b ) , (A4) 

A lb = —± ^Pa-^b + ^a-^b) , (A5) 
Kb = l^a-^b-^ b -A), (A6) 

and {1, ei, e%} are the basis for the field of quaternions |||, which can be represented by two by two matrices, 

1=(J?), e 1= (_°. -;), e^^- 1 ), e 8 =(-' J). (AT) 

Every n x n quaternion matrix A can be represented by a 2n x 2n complex matrix C(A), and the following relations 
hold: det 2 A = detC(A) and 2TiA = TtC(A), where for the quaternion matrix the trace selects the real part 
(component of 1) at the end. Using these properties, the quartic terms in the fermions can then be written as 

T^V'a ■i'bi'b-^a- J^^a ' ^b ' V>a = -^TrA 2 , (A8) 

where the trace also selects the scalar part of the quaternion. These quartic terms can be decoupled by a Hubbard- 
Stratonovich transformation, 

cxp ~ 2a^™ 2 = / DQ exp (~t Tr q2 + iTT A Q ) ' (A9) 

The matrix Q is again a quaternion n x n matrix, 

3 

Q ab = Q° b l+]TQ: h e s , (A10) 

3=1 

with Q° real and symmetric and Q s real and antisymmetric, so that Q has 2n 2 — n degrees of freedom. Such matrices 
are called self-dual real quaternion. The fermions can now be integrated out, and we obtain Eqs.(p7j)-(p8|). 

APPENDIX B: 

In this Appendix we calculate the contributions to the index distribution of saddle point (SP) solutions different 



fromEq. (|32 



for energies in region I. 



A general SP solution reads, 
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Q sp = diag 



(«) 



» j«) 



, z_ 



JO JO 



.,4° 



(Bl) 



V 



"+-p+ 



71— —p- 



For p± ^ Q,n± these SP are not invariant under the unbroken replica symmetry Sp(n + ) x Spin-) of Eq. (p7|), and 
have therefore been called replica symmetry broken (RSB) solutions UM, even though the symmetry is subsequently 
restored by zero modes. The action (p8h at the saddle-point solution (|Bl|) reads, after taking the replica limit, 

iNa(E) [ft- 2ir{p + - p-)] . 

Lets consider now the fluctuations Jl5j . There are sixteen different normal modes. They are labeled by the pair of 
indices (a, a), where the index a = ±1 indicates the upper and lower block, and the index a = ±1 indicates the 
sub-block with solution z a . The eigenvalues are thus denoted by u>r a ,a)(a',tT')- There are three kinds of fluctuations: 

(i) Massive modes: these correspond to a — a' for any a and a' , with eigenvalues wt a <r)(ao J ) = (1 + 

(ii) Zero modes: for a = —a' and a = a' . They are present for any RSB solution. 

(iii) Soft modes: for a = —a' and a = —a 1 , with eigenvalue <jJi a o}i— a -a) = a & e / 7: po(E). 

Lets consider SP's with p + = P— = p, which respect the symmetry of the problem under simultaneous complex 
conjugation and inversion of /i. For these, after integrating out the fluctuations, we obtain 



lim Vl + Vl_ [irp 2 (E)] p G(»,E) , 

%±—>±fj,/2-rr T 



(B2) 



where G(fi,E) Q is the result from the RS solution ( |37| ) (with e ~ 1/N 1 s , see Eq.(^8|)) and the volume of the zero 
mode manifolds V£ ± is given by 



vi=[^pl{E)] 



I 2p(n-p) 



F p 



T(l + n) 



rim 



r(i + 2j) 



T(l+p)T(l + n-p) fJL r[l + 2(n-j + l)]' 



(B3) 



We need now to determine the zero- modes volume in the replica limit. Using the property of the Gamma function 
r(z)r(l — z) — 7r/sin7T2:, and noting that n± — > ±/x/2tt ~ N/A(E), so we want the large \n±\ limit, we find that 



F p 



(-1)P 2 7T- 4 P 2 ^ 



r 2 (i+ P ) 



nr 2 (l + 2i)l (a^^sin^ 



3=1 



2tt/ 



(B4) 



We can now use Eq. (]2(j) to obtain the contribution to the index distribution of the RSB solutions. For the simplest 
one p = 1 we get, 



P i (K,E)=[An 2 p 2 (E)] 



V i"A 3 (£;) 

[K - Na(E) - 1] exp [ 



[K - Na(E) + 1] exp - 



[K - Na(E) - 1] 



[K - Na(E) + if 
2A(E) 



2A(E) 



(B5) 



For \E\ < 2 this is an 0(1/ log N) contribution to the distribution ( |39| ) obtained from the RS solution. In contrast to 
the case of correlation functions Jl3j ], contributions from higher RSB saddle-points do not vanish, but give contributions 
decreasing by powers of logiV. 

Lets turn now to the external regions \E\ > 2. The RSB solutions here are those with p + — (n_ — p_) = p > 0. 
Evaluating ( |28| ) in these SP's we find that their contributions are suppressed by 



exp 



-pN 



\E\\fW 



log 



£2 

2 



\E\ 



as was already found in |l( 
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